# Make temperature plot

source("code/01-importdata.R")

# Load global and US temperatures

global_temp <- read_csv("data/global_temp.csv") %>% 
  filter(Year > 1899,
         Source == "GISTEMP") %>% 
  rename(year = Year, temp = Mean)

udel_temp <- read_csv("data/udel_temp.csv") %>% 
  select(iso3c = iso, year, temp = UDel_tmp) %>% 
  filter(year>1899,
         iso3c != "ATA") %>% 
  drop_na() %>% 
  arrange(iso3c, year) # %>% mutate(demeaned = resid(lm(temp ~ as.factor(iso3c), data = .)))

# Define deviation from common mid century mean
temp_baseline <- udel_temp %>% 
  group_by(iso3c) %>% 
  filter(year > 1950, 
         year < 1981) %>% 
  summarise(mean_mid_century = mean(temp, na.rm = T)) %>% 
  ungroup()

temp_df <- full_join(udel_temp, temp_baseline, by = "iso3c") %>% 
  mutate(demeaned = temp - mean_mid_century)

## Standard deviations
# Individual country standard deviations, and the standard deviation of those
temp_df %>% 
  filter(year<1980) %>% 
  group_by(iso3c) %>% 
  summarise(sd = sd(temp, na.rm = T)) %$% 
  mean(sd, na.rm=T) # 0.4828489

global_temp %>%
  filter(year<1980) %>% 
  mutate(sd = sd(temp)) %$% 
  mean(sd, na.rm=T) # 0.169


## Sample countries for background tracing
# set.seed(2021)
# sample_set <- sample(x = temp_df$iso3c, size = 183, replace = F) # was 10
# 
# temp_bg <- temp_df %>% 
#   filter(iso3c %in% sample_set)

## Plot

# Parameters
text_size <- 8
background_countries <- "grey70"
panel_lines <- "grey60"

# Plot figure 1
p_temperature_sigma_oneway <- temp_df %>% 
  ggplot(mapping = aes(x = year, y = demeaned)) +
  labs(x = "", y = expression(paste("Deviation (", degree, "C) from average annual local temperature", sep = ""))) +
  geom_rect(aes(xmin=1990, xmax=2017, ymin=-10, ymax=10), alpha = 0.1, fill = "grey90") +
  scale_x_continuous(breaks = seq(1900,2020,20), limits = c(1900, 2030)) +
  scale_y_continuous(breaks = seq(-2,3,1)) +
  coord_cartesian(xlim=c(1900,2020), ylim = c(-2.35, 3.4)) +
  theme(legend.position = "none",
        panel.background = element_rect(fill=NA), 
        panel.grid.major.x = NULL,
        panel.grid.major.y = NULL, # element_line(colour = "grey80"),
        panel.ontop = FALSE,
        text=element_text(size=18),
        axis.text = element_text(size = rel(1.25)),
        axis.line = element_line(size = 1, colour = "black")) +
  # Background countries
  geom_line(aes(group = iso3c),
            size = 0.25, color = background_countries) +
  geom_segment(x = 1900, xend = 2017, y = 0, yend = 0, color = "grey40") +
  geom_segment(x = 1900, xend = 2017, y = 0.48, yend = 0.48, color = panel_lines) +
  geom_segment(x = 1900, xend = 2017, y = 2*0.48, yend = 2*0.48, color = panel_lines) +
  geom_segment(x = 1900, xend = 2017, y = 3*0.48, yend = 3*0.48, color = panel_lines) +
  geom_segment(x = 1900, xend = 2017, y = 4*0.48, yend = 4*0.48, color = panel_lines) +
  geom_segment(x = 1900, xend = 2017, y = 5*0.48, yend = 5*0.48, color = panel_lines) +
  geom_segment(x = 1900, xend = 2017, y = 6*0.48, yend = 6*0.48, color = panel_lines) +
  geom_segment(x = 1900, xend = 2017, y = -0.48, yend = -0.48, color = panel_lines) +
  geom_segment(x = 1900, xend = 2017, y = -2*0.48, yend = -2*0.48, color = panel_lines) +
  annotate("text", label = "Observation period", x = (2017+1990)/2, y = 3.4, vjust=0, size = text_size) +
  # Global trend
  annotate("text", label = expression(paste("1", sigma)), x = 2018, y = 0.48, size = text_size, hjust = 0, color = panel_lines) + # 0.48 is the country-level sd
  annotate("text", label = expression(paste("2", sigma)), x = 2018, y = 2*0.48, size = text_size, hjust = 0, color = panel_lines) +
  annotate("text", label = expression(paste("3", sigma)), x = 2018, y = 3*0.48, size = text_size, hjust = 0, color = panel_lines) +
  annotate("text", label = expression(paste("4", sigma)), x = 2018, y = 4*0.48, size = text_size, hjust = 0, color = panel_lines) +
  annotate("text", label = expression(paste("5", sigma)), x = 2018, y = 5*0.48, size = text_size, hjust = 0, color = panel_lines) +
  annotate("text", label = expression(paste("6", sigma)), x = 2018, y = 6*0.48, size = text_size, hjust = 0, color = panel_lines) +
  annotate("text", label = "Mean", x = 2018, y = 0*0.48, size = text_size, hjust = 0, color = panel_lines) +
  annotate("text", label = expression(paste("-1", sigma)), x = 2018, y = -0.48, size = text_size, hjust = 0, color = panel_lines) +
  annotate("text", label = expression(paste("-2", sigma)), x = 2018, y = -2*0.48, size = text_size, hjust = 0, color = panel_lines) +
  # Global observed
  geom_line(data = global_temp, 
            aes(x = year, y = temp),
            size = 1.5, color = "black") +
  geom_point(data = global_temp, 
             aes(x = year, y = temp),
             size = 2, color = "black") 
ggsave(plot = p_temperature_sigma_oneway, "text/images/temperature_sigma_oneway.png", units = "in", height = 8, width = 14, dpi = 320)



# Plot figure 1B: twoway

p_temperature_sigma_twoway <- udel_temp %>% 
  mutate(resid = resid(lm(temp ~ as.factor(iso3c) + as.factor(year), data = .))) %>% 
  ggplot(., aes(as.numeric(year), as.numeric(resid))) +
  # geom_line() +
  labs(x = "", y = expression(paste("Deviation (", degree, "C) from national and global temperature", sep = ""))) +
  geom_rect(aes(xmin=1990, xmax=2017, ymin=-10, ymax=10), alpha = 0.1, fill = "grey90") +
  scale_x_continuous(breaks = seq(1900,2020,20), limits = c(1900, 2030)) +
  scale_y_continuous(breaks = seq(-2,3,1)) +
  coord_cartesian(xlim=c(1900,2020), ylim = c(-2.35, 3.4)) +
  theme(legend.position = "none",
        panel.background = element_rect(fill=NA), 
        panel.grid.major.x = NULL,
        panel.grid.major.y = NULL, # element_line(colour = "grey80"),
        panel.ontop = FALSE,
        text=element_text(size=18),
        axis.text = element_text(size = rel(1.25)),
        axis.line = element_line(size = 1, colour = "black")) +
  geom_line(aes(group = iso3c),
            size = 0.25, color = background_countries) +
  annotate("text", label = "Observation period", x = (2017+1990)/2, y = 3.4, vjust=0, size = text_size) +
  # Global observed
  geom_line(data = global_temp, 
            aes(x = year, y = temp),
            size = 1.5, color = "black") +
  geom_point(data = global_temp, 
             aes(x = year, y = temp),
             size = 2, color = "black") 
ggsave(plot = p_temperature_sigma_twoway, "text/images/temperature_sigma_twoway.png", units = "in", height = 8, width = 14, dpi = 320)
